# 通过人类相互作用组将多酚靶点与致病蛋白联系起来
import pandas as pd
import networkx as nx
# from matplotlib import pyplot as plt 
#数据保存
# data_file =r'./EECG.csv'
# data_file=r'./'+disease+'.csv'
#数据保存到csv文件，使用逗号分割
def save_data(filename,data_to_write):
  with open(filename,"a") as f:
    f.write(data_to_write)

# 读入多酚名称与靶点蛋白
polyphenol = pd.read_csv('./data/PolyphenolProteinInteractions.csv', 
                         dtype={'entrez_id': int})
#   chemical  pubchem_compound_id  entrez_id  symbol
# 0   butein            5281222.0       6718  AKR1D1

#多酚名称 通过set进行 去重 并按照名称排序
#set  1：去重
#     2: 排序
polyphenol_list = set(polyphenol['chemical'])

# 读入疾病分类与节点蛋白
dg = pd.read_csv('./data/GenesDisease.csv', dtype={'entrez_id': int})

# 节点蛋白使用set 排序并 去重
dg_list = set(dg['disease'])

#            disease  entrez_id
# 0  kidney diseases      79663

# 读入人类相互作用组，并构成网络
dt = pd.read_csv('./data/HumanInteractome_v2017.csv',
                 dtype={'EntrezA': int, 'EntrezB': int})
#    EntrezA  EntrezB
# 0        1      310

#根据人类相互作用组构建网络
#使用networkx的函数，根据前面读取到的 人类相互作用组 作为 “边” 信息，构成网络
# A -->  B 
# B -->  C
# A -->B --> C   ===>  A --> C 
G = nx.from_pandas_edgelist(dt, 'EntrezA', 'EntrezB')

# 只考虑相互作用组中最大连通分量(LCC)
largest_cc = max(nx.connected_components(G), key=len)
g = G.subgraph(largest_cc)


# nx.draw_networkx(g)
# plt.show()
#通过封装成函数
# 参数 是 多酚 、疾病蛋白 ，子图
# 作用：根据 多酚 、疾病蛋白分别 与 子图 的交集 计算 dc
def caculate_dc(chemical,disease,g_sub):
    # 多酚到疾病的距离
    # chemical = '(-)-epigallocatechin 3-o-gallate'
    # 多酚名称与靶点蛋白
    # line 61  67 根据 前面的表格 找出药物名对应的 id 值 ，因为图通过ID构建 
    target_nodes = polyphenol[polyphenol['chemical'] == chemical]['entrez_id']
    # print('{}与{}个蛋白质有关'.format(chemical, len(target_nodes)))

    # disease ='glucose metabolism disorders'
    # disease ='nervous system diseases'

    disease_nodes = dg[dg['disease'] == disease]['entrez_id']
    # print('{}与{}个蛋白质有关'.format(disease, len(disease_nodes)))

    # 两组蛋白间可通过其它组的蛋白联接
    # g_sub = g
    #取交集
    t_nodes = set(target_nodes) & set(g_sub)   #a
    d_nodes = set(disease_nodes) & set(g_sub)  #b

    min_lengths = []
    mean_lengths = []
    #遍历源点到目标点集合
    #将最短路径 追加 到  t_length 中，找出 min  mean 
    for node_from in t_nodes: #a 1 2...
        t_length = []
        for node_to in d_nodes: #b  1 2 3 ...
            if (nx.has_path(g_sub, node_from, node_to)):
                t_length.append(nx.shortest_path_length(g_sub, node_from, node_to))# 最短路径长度

        t_length = pd.Series(t_length)
        min_lengths.append(t_length.min())
        mean_lengths.append(t_length.mean())

    min_lengths = pd.Series(min_lengths)
    mean_lengths = pd.Series(mean_lengths)
    #返回计算结果
    return min_lengths.mean(),mean_lengths.mean()
# 保留小数点后三位有效数字输出   line 116
# print('从{}到{}的d_min为： {:.3f}'.format(chemical, disease,
#                                           min_lengths.mean()))
# print('从{}到{}的d_mean为： {:.3f}'.format(chemical, disease,
#                                            mean_lengths.mean()))





def main():
    #选定一种多酚
    chemical='(-)-epigallocatechin 3-o-gallate'
    data_file=r'./out/'+chemical+'.csv'
    #已经计算完成的计数
    count = 0
    #表头
    save_data(data_file,'chemical'+','+'disease'+','+'min_len'+','+'mean_len\n')
    #遍历 疾病蛋白 集合
    for disease in dg_list:
        count+=1
        [min_len,mean_len]=caculate_dc(chemical,disease,g)
        # print(chemical+'->'+disease +':',min_len,mean_len)
        print('[{}/{}]  {}    :   {}    d_min：{:.3f}  ,d_mean:{:.3f}  '.format(count,len(dg_list),chemical, disease,min_len,mean_len))
        save_data(data_file,chemical+',"'+disease+'",'+str(min_len)+','+str(mean_len)+'\n')

#保证从main函数开始运行
if __name__ == '__main__':
    main()